Sheared force- networks: anisotropies, yielding and geometry 
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A scenario for yielding of granular matter is presented by considering the ensemble of force 
networks for a given contact network and applied shear stress r. As r is increased, the probability 
distribution of contact forces becomes highly anisotropic, the difference between average contact 
forces along minor and major axis grows, and the allowed networks span a shrinking subspace of 
all force-networks. Eventually, contacts start to break, and at the yielding shear stress, the packing 
becomes effectively isostatic. The size of the allowed subspace exhibits simple scaling properties, 
which lead to a prediction of the yield stress for packings of arbitrary contact number. 
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Granular media, foams and emulsions are amorphous 
materials which can jam and then sustain a certain 
amount of shear stress before yielding 0, S 0> IS- 
If one slowly increases the applied shear stress and fol- 
lows the evolution of contact forces and grain locations 
for such systems, one encounters a rather complex set 
of phenomena. Firstly, before the system yields as a 
whole, there are non-adiabatic precursor events such as 
local rearrangements and microslip 0, Qi ■ Secondly, 
the interparticle contact forces in these systems are or- 
ganized into high ly heterogeneous and fragile force net- 
works IS [H O El El - see Fig. n Unraveling these 
shear-induced phenomena and their impact on macro- 
scopic unjamming has remained a great challenge. 

While the contact forces evolve to satisfy the applied 
stresses, the selection of a specific force network for a sin- 
gle numerical experiment hinges on microscopic details 
and packing history. On the other hand, features like con- 
tact force probabilities have turned out to be relati vely 
insensitive to these subtle details [H E [Ji', li, Ig, 
which suggests a purely statistical approach. In this Let- 
ter we characterize granular packs under shear stress, by 
studying ensembles of force networks for fixed contact 
networks JJe J^, JJ„, 20, ,21j. This approach is based 
on the fact that in jammed systems there are more con- 
tact forces than force balance equations - the ensemble 
simply consists of all those force networks for which the 
contact forces are repulsive, balance on every grain and 
satisfy global stress constraints, while keeping the geom- 
etry fixed. This provides a novel access to the statisics 
of force networks under shear stress, in which the roles 
of fabric and force anisotropy are separated explicitly. 

We consider ensembles of sheared force networks for 
frictionless disks in two dimensions, for contact numbers 
z ranging from the lower limit z « Zc = 4 (isostatic [2^ ) 
to z = 6 (strongly hyperstatic). We recover a number of 





FIG. 1: (a) Overview of the z-t parameter space, where z and 
T denote coordination number and shear stress respectively. 
For a given z, force networks cease to exist beyond a maximum 
Tm (dashed line). The dots indicate the values of r and 2 in 
panels (b-d). (b-d) Parts of the force networks, where line 
thicknesses represent the strength of the contact forces. These 
networks become highly anisotropic under shear stress. 



experimental features, in particular a transition to yield- 
ing: force networks cease to exist beyond a critical yield 
stress T„j. The boundary of the grey zone in Fig. la in- 
dicates that this maximum yield stress strongly depends 
on the coordination number. 

In this paper we characterize the physics of sheared 
force networks at three levels of detail, (i) P(j,{f), the 
contact-angle (0) resolved probability density for the con- 
tact forces /, is a natural extension of the overall force 
distribution P(/) [H E El H [13 , and acts as a sen- 
sitive probe for anisotropy due to shear, (ii) The an- 
gle resolved average force /(0) is possibly the simplest 
characterization of shear-induced anisotropy. We find 
that it evolves sinusoidally as /(0) = 1 -I- 2t sin(20) for 
small T, but that higher harmonics develop for larger 
shear stresses. Bounds on r™ are obtained by requiring 
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FIG. 2: (a) Illustration of our geometry, showing the direc- 
tion of shear stress (arrows), contact angle and major and 
minor axes, (b) Angle resolved force distributions for z = 5. 
P<l>{f) becomes increasingly modulated with contact angle (p 
for increasing shear stress r. (c) P^{f) varies qualitatively 
with cj} (here r = 0.2). 
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FIG. 3: Average contact force as function of contact angle 
4> and applied shear stress r. (a) f{cf)), for z = 5 develops a 
sinusoidal modulation when r is increased. /(</>) can be fitted 
well by an expression of the form 1 + 2Tsin(0) — 62 cos(4(;/>); 
62 slowly increases with r as shown in panel c. (b) /(<;/>), 
for z = & develops a strong second harmonic for large r and 
approaches the limiting curve given by Eq. Q (dashed curve), 
(c-d) Corresponding dependence of &2 on r. 



that /((/)) should not become negative, (in) The yield- 
ing curve of Fig. la is due to a vanishing volume of the 
high-dimensional space of allowed networks F,- . Charac- 
terizing the size of Ft- by the average distance L between 
random pairs of force networks, we find that L and r are 
approximately related as {L/Lmf + [T/TmY = 1- Both 
Lm and exhibit simple scaling laws with contact num- 
ber z, z — Zc and particle number N , resulting in a general 
prediction for the maximal shear stress that a given 
granular packing can sustain before yielding. 

Force ensemble - The numerical results have been 
obtained by a recently developed ensemble technique 
mill III 113, El- The input of the ensemble consists 
of a fixed contact geometry of a 2D packing of = 1024 
frictionless disks of radii Ri with centers and coordina- 
tion number 2: > Zc in a volume V , which we generated 
from molecular dynamics simulations of a 50:50 binary 
mixture of particles with size ratio 1.4 that have a purely 
repulsive Lennard-Jones interaction. Different densities 
were used to obtain different z. The system is not sheared 
and the resulting contact networks are isotropic H ^3 ■ 
These contact networks are then kept fixed, and the pos- 
itive interparticle forces between particles i and j, fij, 
are seen as degrees of freedom which satisfy mechanical 
equilibrium, restricted by the macroscopic stress aa/3'- 

(1) 

In this picture there are {zN/2) degrees of freedom (con- 
tact forces) constrained by (2iV -I- 3) equations, leading, 
for z > 4, to an ensemble of force networks that form a 
high-dimensional convex subspace Ft- [rRll8lIl9| . In this 
paper we choose the coordinates and pressure such that 



o'xx — <^yy = 1/2, and consider the dimensionless shear 
stress T — Oxyl <^xx (equivalent to the relative deviatoric 
stress) (see Fig. 2a). 

Angle resolved force distributions ~ To characterize 
the anisotropic force networks, we introduce here the 
contact-angle resolved force distribution P^pif)- This dis- 
tribution is normalized such that / dfP^{f) = 1 for all 
0, and I/tt J dcf) J df[fP^{f)] = 1. Note that, in gen- 
eral, the average force for contacts in a certain direction, 
f{(l)) = J df [fP^if)], is not equal to one. 

Figure 2b illustrates that for sheared systems P^pif) 
modulates with (j) - this modulation has its extrema along 
major and minor axes. Many contact forces along the 
minor axis {(j) ^ 37r/4) evolve towards zero for increasing 
shear, effectively breaking these contacts 17j and lead- 
ing to a (S-peak at / = (black area in Fig. 2b4). Fig. 2c 
shows examples of Pcf,{f) and illustrates that for sheared 
systems, the qualitative shape of the force distribution 
varies with contact angle. In particular for large stresses, 
Pct>{f) becomes completely monotonic along the minor 
axis but remains peaked along the major axis. Earlier 
experimental and numerical findings that the total force 
distribution, averaged over all orientations, evolves from 
a peaked to a monotonically decreasing function with in- 
creasing r are thus seen to result from averaging over 
strongly different ^^(7) 011 El 13111 ■ 

Analytical bounds on t - The most basic manifesta- 
tion of stress anisotropy, however, is the modulation of 
the average force, f{4>), as a function of the contact ori- 
entation. This effect is clearly visible in FigSjL and has 
only recently been accessed experimentally [liJ, |l3| ■ In 
Fig. 3 we therefore show examples of f{(t>) for various 
stresses and contact numbers, as obtained by the ensem- 
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ble. For the strongly hyperstatic case (Fig. 3b, d), it is 
clear that the anisotropy is limited by the requirement 
that f{(j)) should definitely remain positive for all (f>. This 
is due to the repulsive nature of the contact forces, which 
requires all fij > 0. We show below that this simple cri- 
terion imposes analytical bounds on the maximum shear 
stress T„i that become increasingly accurate for strongly 
hyperstatic packings. 

This bound can be computed from the probabilistic 
version of Eq. ^ for isotropic contacts, 

= "jT^ j d(t) f{(j))nanp , (2) 

where r is the average particle radius '5^, while 
(rix ,ny) = (cos (f>, sin (j)) ; in the remainder we set the pref- 
actor fNc/V equal to unity. We Fourier expand /((/>) as 
^ Uk sin 2k(j) + J2 cos 2/c(/), and anticipate that only 
odd k sine terms and even k cosine terms are compatible 
with the symmetry of a simple shear 0. Eq. © then 
yields 

/((/>) = 1 + 2t sin 2<?!) - 62 cos 4<?!) +••• , (3) 

where the coefficients of the higher order terms 62 , • ■ ■ 
are independent of the stress tensor (Fourier modes with 
k > 2 yield zero upon integration). The role of the higher 
order terms is limited 24]: our numerics show a signifi- 
cant contribution to /(</») for large r and z for k — 2 only 
(Fig. 3b, d). Truncating Eq. ||2J) at second order yields 
that an upper bound is reached for T = Tm = 2/3: 

fmax{(f>) = l-H4/3sin20 - 1/3 COS 40 (4) 

By extending the expression ^ to include fabric 
anisotropy and frictional forces, and following a similar 
line of reasoning, bounds can also be obtained for the 
completely general case |2^ . 

Fig. 3b, d illustrate the relevance of this bound for 
strongly hyperstatic packings: for z — 6 the maximal 
stress is close to 2/3, while /(0) approaches the limit- 
ing form Eq. Q), indicated by the dashed Hne. Closer to 
isostaticity, however, the system has much less freedom to 
accommodate to the shear stress, and Tm is much smaller 
- as shown in Fig. 3a, c, / then does not approach zero 
yet. In the isostatic limit there are no adjustable degrees 
of freedom left, the force space F-r shrinks to a single 
point, and tends to zero. In general, the question is 
thus how Tjn depends on z (Fig. la). 

High- dimensional ensemble - The yielding curve 
Tm{z) can be determined best from the geometric proper- 
ties of this high-dimensional force space. We have quan- 
tified the " size" of F,- by the Euclidian distance between 
random realizations (force networks), for given z and r 
^jHO]. While the distances between random points in a 
low dimensional space are broadly distributed, this dis- 
tribution becomes increasingly sharply peaked for higher 




-0.7 0.0 0.7 




FIG. 4: (a) The linear size L of the force space Ft as a func- 
tion of shear stress r for 2 as labeled (dots) are well fitted by 
a relation of the form [L/Lm.)^ + (j jTrnf" = 1 (solid curves), 
(b) Sketch of how force space F (represented schematically 
by a sphere) intersected by the constant-r hyperplanes yields 
the force spaces and Fra with sizes L\ and L^. (c) L„i 
obtained by fitting Eq. Q (circles) , compared to ^ N{z — Zc) 
(line), (d) Maximal shear stress Tm obtained by fitting Eq. Q 
(circles) is well approximated by 2{z — Zc)/z (line), (e) Zctt 
drops sharply when approaching Tm- 

dimensional objects as is the case here. The average dis- 
tance L defined via 

L\z,T)^{E,,if,,-n^)') , (5) 

thus serves as an effective measure for the size of F,-. 
Here the brackets denote an average over the random 
pairs of force networks {fij} and {fij} 

In Fig. 4 we show our main findings for the main prop- 
erties of the force space with r. First, L(t, z) can, sur- 
prisingly, be fitted by a simple relation of the form 

{L/L,nf + (r/r„02 ^ 1 ^ (5) 

which becomes particularly accurate for z < 5 (Fig. 4a). 
In addition, as shown in Fig. 4c-d, L„i and Tm approach 
simple scaling laws: 

r„«2i^, (7) 

z 

Lm - VN{z - z,) . (8) 

The scaling relations Eqs. H6I8|I can be interpreted geo- 
metrically, keeping in mind that high-dimensional objects 
can be quite counter-intuitive (Fig. 4b). Let us consider 
the ensemble F that is obtained by applying all force bal- 
ance equations and Oxx—cfyy = 1/2, but leaving a^y, and 
thus T, undetermined. F is a convex body of dimension 
D — {z~ Zc)N /2 — 2 that has zN /2 facets; each facet cor- 
responds to a certain contact force being zero. One can 
obtain F,- from intersection of F with the codimension- 
one hyperplane given by the linear constraint axy = t/2 
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(Eq. Due to symmetry, the most "central" intersec- 
tion is obtained for r = 0, while for larger values of t, 
the intersection is less centered and F,- is smaller. 

In fact, Eq. © can be recovered analytically. The 
scaling Lm esc ^/Ij with the dimension D = N/2{z — Zc) 
of Ft is a common feature of high-dimensional convex 
spaces |27f . If we consider a simple approximation for F^. 
by ignoring the force balance equations but only requiring 
that all fij are positive and have (fij) = 1, we obtain 
that = 2D = N{z - Zc) Note that the relation 

Tm ~ 2(z — Zc)/z, implies that is proportional to the 
ratio between the dimension of F^ and the number of its 
facets zN/2. We have not been able to come up with a 
convincing argument for this scaling, however. 

Finally, as t approaches its maximum r^, the space 
Ft- shrinks and L ^ 0, so that at T = Tm, consists of a 
single point. The effective contact number ZeS, which is 
defined by considering contacts broken when their force 
drops below a fixed small threshold ^2^, stays constant 
over most of the range of r, but sharply drops to 
approaches r™ (Fig. 4e). 

Outlook - We have proposed an ensemble theory for 
force networks under shear, by exploring solutions of me- 
chanical equilibrium for packings of fixed contact geom- 
etry. The ensemble approach provides a great concep- 
tual simplification with respect to full numerical simu- 
lations, as it steps aside the intricate evolution of the 
contact network. Yet, it captures recently measured sta- 
tistical properties, such as f{(f>) and the evolution of P(/) 
Furthermore, it provides an alternative descrip- 
tion of yielding phenomena, in terms of a vanishing vol- 
ume of the force phase space. Consistent with existing 
numerical simulations we found that the maximum 

shear stress t^, strongly depends on the coordination 
number of the packing. This dependence can be under- 
stood in terms of the geometry of the force space, and 
obeys a simple scaling law Eq. (0. The ensemble thus 
provides a new perspective for soil mechanics, in which 
relations between the macroscopic effective friction and 
micromechanical properties (density, coordination num- 
ber, texture,... ) play a central role. More generally, it 
suggests a route along which the unjamming by shear of 
a broad range of disordered media may be understood. 
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